Load data from files:

library("tidyverse")
library("cowplot")
library("readxl")
library("openxlsx")
library("plotly")
library("rcartocolor")
library("tidytext")


### Load Kynnetec AgroTrak pesticide data (proprietary):
AgroTrak.CornSoy.dat <- read_excel("AgroTrak_CornSoy_v1.xlsx",
                          skip = 5) %>%
  select(Year, Crop,
         Type = "Pesticide Type",
         Active.Ingredient = "Active Ingredient",
         Volume.lbs = "AI volume\r\n(lb)",
         Base.Area.Treated = "Base Area Treated\r\n(acres)") %>%
  mutate(ai = str_replace_all(Active.Ingredient, " \\(.\\)", "")) 
glimpse(AgroTrak.CornSoy.dat)

### Load concordance file:
srcTab <- read_csv("cclistInEcotox.csv")
glimpse(srcTab)

### Load honey bee toxicity values from ecotox export:
ecotox.dat0 <- 
  read_excel("ECOTOX-Terrestrial-Export_20211207_plusManualAdd.xlsx",
             sheet = "Terrestrial-Export") %>%
  rename(cas.number = "CAS Number", chemName = "Chemical Name") 
### Pare down ecotox values to those likely to be used in this analysis
### (leaving the full ecotox.dat0 just in case needed later),
### and join with the concordance list containing AgroTrak names:
ecotox.dat <- ecotox.dat0 %>%
  select(source,
         cas.number, 
         ConcType = `Conc 1 Type (Author)`,
         Effect, 
         LifeStage = `Organism Lifestage`,
         ExposureType = `Exposure Type`, 
         Endpoint, 
         ObservedDuration.d = `Observed Duration (Days)`,
         ObsRespMean = `Observed Response Mean`,
         ToxUnit = `Observed Response Units`) %>%
  filter(Endpoint %in% c("LD50", "NOEL")) %>%
  right_join(srcTab) %>%
  select(-inEcoTox, -inManualAdd, -CompTox.PREFERRED_NAME)
glimpse(ecotox.dat)

Create a ranked list of pesticides by volume applied in AgroTrak:

AgroTrak.Ranked <- AgroTrak.CornSoy.dat %>%
  group_by(Type, ai) %>%
  summarize(Volume.lbs = round(sum(Volume.lbs, na.rm = TRUE)),
            Base.Area.Treated.sum = round(sum(Base.Area.Treated, 
                                              na.rm = TRUE)),
            nYears = length(unique(Year)),
            Base.Area.perYear = round(Base.Area.Treated.sum / nYears),
            Volume.perYear = round(Volume.lbs / nYears)) %>%
  arrange(Type, -Volume.lbs)

Adult Contact LD50

AdultContactLD50 <- ecotox.dat %>%
  filter(LifeStage == "Adult" &
           ExposureType == "Dermal" & 
           Endpoint == "LD50") %>%
  group_by(ai) %>%
  mutate(checkAI = any(grepl("AI|ai|ae", ToxUnit)),
         hasAI = grepl("AI|ai|ae", ToxUnit)) %>%
  group_by(ai, hasAI) %>%
  mutate(min.d.grp = min(abs(2 - ObservedDuration.d) + 2, 
                         na.rm = TRUE)) %>%
  filter(case_when(checkAI == TRUE ~ 
                     hasAI == TRUE & ObservedDuration.d <= min.d.grp + 2,
                   checkAI == FALSE ~
                     ObservedDuration.d <= min.d.grp + 2)) %>%
  ungroup() %>%
  select(-checkAI, -hasAI, -min.d.grp) %>%
  group_by(ai) %>%
  slice(which.min(ObsRespMean))
glimpse(AdultContactLD50)
unique(AdultContactLD50$ConcType)
unique(AdultContactLD50$ToxUnit) 
AgroTrak.AdultContactLD50.Rank <- left_join(AgroTrak.Ranked, AdultContactLD50)
#write.xlsx(AgroTrak.AdultContactLD50.Rank,
#          file = "Draft_AdultContactLD50.xlsx", asTable = TRUE)

AdultContactLD50.summary <- AgroTrak.AdultContactLD50.Rank %>%
  group_by(Type) %>%
  mutate(totalVol = sum(Volume.lbs),
         pctVol = Volume.lbs / totalVol * 100,
         toxVal = !is.na(ObsRespMean)) %>%
  #select(Type, ai, Volume.lbs, totalVol, pctVol, toxVal)
  group_by(Type, toxVal) %>%
  summarize(toxPct = sum(pctVol)) %>%
  pivot_wider(names_from = toxVal, 
              values_from = toxPct) %>%
  rename("missing Tox value" = `FALSE`,
         "Tox value present" = `TRUE`)
knitr::kable(AdultContactLD50.summary, digits = 1)
Type missing Tox value Tox value present
Fungicide 6.0 94.0
Herbicide 3.4 96.6
Insecticide 8.3 91.7

Analysis of Adult Contact LD50 RQs for Corn & Soybean

List of active ingredients that made up greater than 1% of a group’s volume in any single year but have no adult dermal LD50 value:

CornSoyAdultContact.dat %>%
  filter(is.na(RQ.pctOfType) & Volume.pctOfType > 1) %>%
  distinct(Type, ai) %>%
  arrange(Type)
## # A tibble: 17 × 2
##    Type        ai               
##    <chr>       <chr>            
##  1 Fungicide   SULFUR           
##  2 Fungicide   TEBUCONAZOLE     
##  3 Fungicide   COPPER SULFATE   
##  4 Fungicide   TETRACONAZOLE    
##  5 Fungicide   FLUSILAZOLE      
##  6 Fungicide   PETROLEUM OIL    
##  7 Fungicide   PHOSPHORIC ACID  
##  8 Herbicide   CYANAZINE        
##  9 Herbicide   ACIFLUORFEN      
## 10 Herbicide   BENTAZONE        
## 11 Herbicide   FOMESAFEN        
## 12 Herbicide   SULFOSATE        
## 13 Herbicide   GLUFOSINATE      
## 14 Herbicide   DIMETHENAMID-P   
## 15 Insecticide TEBUPIRIMPHOS    
## 16 Insecticide BIFENTHRIN       
## 17 Insecticide ZETA-CYPERMETHRIN

Figures

Total pesticide volume applied.

Pesticide volume applied per treated acre.

Summed honey bee risk quotient values for all pesticides.

Summed honey bee risk quotient values per treated acre.

Same RQ plots, faceted & scaled to see fungicide & herbicide trends.

## `summarise()` has grouped output by 'Crop', 'Type'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'Year', 'Crop', 'Type'. You can override
## using the `.groups` argument.
## Rows: 1,248
## Columns: 7
## Groups: Year, Crop, Type [116]
## $ Year           <dbl> 1998, 1998, 1998, 1998, 1998, 1998, 1998, 1998, 1998, 1…
## $ Crop           <chr> "Corn", "Corn", "Corn", "Corn", "Corn", "Corn", "Corn",…
## $ Type           <chr> "Herbicide", "Herbicide", "Herbicide", "Herbicide", "He…
## $ ai.top         <chr> "2,4-D", "ACETOCHLOR", "ALACHLOR", "ATRAZINE", "DICAMBA…
## $ BeeREX.RQ      <dbl> 446199.10, 463729.70, 591170.16, 1927358.29, 165447.57,…
## $ BeeREX.RQ.acre <dbl> 0.07176011, 0.02434560, 0.15482690, 0.03452120, 0.01014…
## $ RQ.pctOfType   <dbl> 7.6780180, 7.9796776, 10.1726227, 33.1652207, 2.8469564…

Soybean Insect Pest Targets:

The “RQ Pest Index” needs some scrutiny…